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Abstract 

Recent studies have revealed much of the mathematical structure of the 
static correlation functions of the XXZ chain. Here we use the results of 
those studies in order to work out explicit examples of short-distance cor- 
relation functions in the infinite chain. We compute two-point functions 
ranging over 2, 3 and 4 lattice sites as functions of the temperature and the 
magnetic field for various anisotropics in the massless regime —1 < A < 1. 
It turns out that the new formulae are numerically efficient and allow us to 
obtain the correlations functions over the full parameter range with arbi- 
trary precision. 
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1 Introduction 



In recent years considerable progress was achieved in the understanding of the static 
correlation functions of local operators in the XXZ spin chain. This concerns mostly 
the ground state correlation functions at zero magnetic field for which a hidden fermi- 
onic structure (outside the free Fermion point A = 0) was identified in [4-6]. The ex- 
ponential of an operator, bilinear in two types of annihilation operators and acting on 
the space of local operators on the spin chain, generates all local correlation functions, 
while the space of local operators on the chain is generated by the corresponding cre- 
ation operators. For the inhomogeneous chain this structure implies the factorization 
of arbitrary correlation functions into sums over products of neighbour correlators. 

A generalization of the exponential form to finite temperatures and finite longi- 
tudinal magnetic field was suggested in [2]. Without magnetic field the same Fermi 
operators as in the ground state can still be used, and the modification of the ground 
state result concerns only certain functions multiplying the annihilation operators in 
the exponent. For finite magnetic field, on the other hand, it seems that the algebraic 
structure must be altered as well and a new Fermi-like operator appears linearly in the 
exponent. 

Here we shall work out the finite temperature formula of [2] in detail and compute 
the longitudinal and transversal spin-spin correlation functions for up to four lattice 
sites. We shall concentrate on the massless regime —1 < A < 1, since the massive 
regime A > 1 requires the calculation of different functions in the analytical part of our 
work, which we postpone to a separate publication. 

2 The density matrix of a chain segment 

We consider the XXZ chain with Hamiltonian 



in the antiferromagnetic regime (J > and A > —1)^. Here the operators a", j = 
—N +1,...,N, a = x, y,z, act locally as Pauli matrices, J is the exchange coupling, 
and A = ch(r|) = \(q + q~ l ) is the anisotropy parameter. 

The XXZ Hamiltonian preserves the z-component of the total spin 




N 



(1) 



j=-N+\ 



N 



QZ _ 1 

Q N - 2 



(2) 



j=-N+\ 



^The results of this and of the following section are valid for all A > — 1 . Later on in sections 4 and 
5 we restrict ourselves to the massless regime — 1 < A < 1 . 
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The equilibrium properties of the chain at temperature T and external magnetic field h 
are determined by the statistical operator 

p N (T,k) = 7777 — 7 s \ /r ■ (3) 

tr^+i^.^e-^-^)/ 7 

We are interested in the system in the thermodynamic limit L = IN — > °°. To 
perform this limit in a sensible way we fix a positive integer n and define 

D n (T,h) = lim tr^+l,. .. -l ,0,n+i,n+2,...,iv pN(T,h) , (4) 

the density matrix of the segment [1 ,n] of the infinite chain. By construction it satisfies 
the reduction property 

tx n D n (T,h)=D n . 1 (T,h). (5) 

Denote by W the space of operators with non-trivial action only on a finite number of 
positive lattice sites. Because of the translational invariance of the Hamiltonian it is 
sufficient to consider this space. For G W we write its restriction to the first n lattice 
sites as 0[i, n ]. Then 

D* TM (0) = lmte w (D„(r,fc)0 [1)B] ) (6) 

is well defined because of (5) and determines the thermal average of the operator in 
a magnetic field h, which we shall also denote (0)r,h- Mathematically (6) defines Dj . 
as a functional on W. Along with D* h we often study its inhomogeous generalization 
(see [2]) which we denote by the same letter. When acting on an operator of length n 
it depends on n so-called inhomogenieties j = 1 , • • • , n. 

In our previous work [2] we conjectured a rather explicit formula for Dj h which 
was inspired in part by [5]. The conjecture says that 

D* Th {Q) =tr(e fil+fi2 (0)), (7) 

where tr is the trace functional defined by 

tr(0) = ...itn ±tr 2 itr 3 ...(0), (8) 

and the Clj : W — > W, j = 1 , 2 are linear operators of the form* 

ai = -^o/ r /r2^2^^ (9a) 
Q 2 = -lini ,f§(p(// i; a)h(^a). (9b) 



^We have slightly changed our notation for operators. Instead of those of [2] we rather employ the 
conventions of the more recent paper [6]. 
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Here fij = ln^-, j = 1,2. As a consequence of taking the limit a^Owe shall only 
need the values of the function Cfl(//i,//2;oc) and of its a-derivative at a = 0. Similarly 
we only need (p(//i;0). All the functions that are necessary for the actual evaluation of 
(9) will be introduced in section 4 below. They depend on the temperature and on the 
magnetic field and describe what we call the physical part of the problem. By way of 
contrast, the operators b(£,cc), c(£,cc) and h(£,oc) are independent of T and h. These 
operators are rational functions of ^ 2 and constitute the algebraic part of the problem. 
We will explain their construction in the next section. The closed contour T encircles 
the poles of the operators at the inhomogeneities ^ (see below) and excludes all other 
singularities of the integrand. 

3 The algebraic part of the construction 

We should first explain that our present status of understanding is different for Q.i and 
^2 in (9). The algebraic part of Qi, built up from the operators b(£,cc) and c(£,cc), 
is the same as for the ground state at vanishing magnetic field. First versions of these 
operators were introduced in [5]. They were optimized and further studied in [4,6]. 
Their analytic and algebraic structure is by now well understood. In particular, they 
satisfy a reduction property and anticommutation relations. The operator h(£,oc) was 
introduced in [2] in order to account for correlation functions odd under spin reversal. 
This operator is less well studied. We plan to further elaborate on its mathematical 
structure elsewhere. Below, when we work out examples of correlation functions, we 
use some of its properties for a — > that we verified by explicit calculations for up to 
four lattice sites. 

In this section we follow to a certain extent the recent paper [6] to which we refer 
the reader for further details. For the definition of the operators b(£, a), c(£, a) we first 
of all generalize the space of local operators. Instead of local operators we consider 
operators of the form 

q 2aS(0) Qj S(0) = I£c^., (10) 

J = — 00 

where is local and a e C. These operators are called quasi-local with tail a. They 
span a vector space W a . The parameter a is introduced in order to regularize the 
various objects introduced below. Some of them are rather singular for a — > 0, and the 
limit is, in general, only well defined for appropriate combinations. 

In [6] the operators b(£,cc), c(£,cc) were defined in two steps. In step one endo- 
morphismsbnyi(£,a) and c a) acting onMjyi = End(Vfc<g>---<E>V/) were defined, 
where every Vj = C 2 is the space of states of a spin in the infinite chain. In step two a 
certain reduction property was used to extend their action inductively to W a . Here we 
shall only review step one, since we will be dealing with small finite chain segments 
anyway, and we refer the reader to [6] for the reduction property. 

The operators bpyj (£, a) and cpyj (£, a) are constructed from weighted traces of the 
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elements of certain monodromy matrices related to t/ 9 (s b). These monodromy ma- 
trices are products of two types of L-matrices with two-dimensional and with infinite 
dimensional auxiliary space, respectively. 

The L-matrices with two-dimensional auxiliary space are directly related to the 
/^-matrix of the six- vertex model, 



R(Q = 



(\ o\ 

(3(0 7(0 

y(Q (3(0 

\0 \) 



(11) 



where 

Fixing an auxiliary space V a isomorphic to C 2 we define L a j(Q — p(QR a j(Q, where 
p(0 is a scalar factor to be specified below. This is the standard L-matrix of the six- 
vertex model. The corresponding monodromy matrix is 

T am (Q=L a 4^ k )...L a ^). (13) 

It acts on V a <8> VJt <8> • • • <8> V/. We are interested in operators acting on M^/j or on 
M a <g).A% n, where M a = End(V a ). Such type of operators are naturally given by the 
adjoint action of operators acting on V* ® • • • <g) V/ or on V a (g) V]t (8) • • • <g> V/. An example 
is the L-operator L a j(Q defined by 

L flJ (0%/]=L flJ (0% ; ^-)(0 (14) 

for j E {&,...,/} and for all e M a ®M^j. These L-operators generate the first 
type of monodromy matrices entering the constructions of b(0 oc), c(0 oc), 

T a (^a)X M = Lat&fe) ■ ■ .K,i(WQ aGZ X[k,i] (15) 

far all X [Jt)/] eM a ®M [k fl. 

A second type of monodromy matrix acts adjointly on an infinite dimensional 
bosonic auxiliary space, more precisely on a representation space of the g-oscillator 
algebra Osc defined in terms of generators a, a*, q ±D and relations 

q D nq- D = q- 1 a , q D H*q- D = qa\ aa* = 1 - q 2D +\ Z*Z=\-q 2D . (16) 

Let Osc a a copy of Osc. We define 
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which acts on Osca <8> Vjt <8> • • • <8> V/. The scalar factor o(Q is a convenient normalization 
and is required to solve the functional equation 



1 



l-^ 2 



(18) 



It also fixes the scalar factor in the definition of L a j(Q, 



p(Q 



For the L-operator (17) we define again its adjoint action 



(19) 



(20) 



for j e {A:, ...,/} and for all X^/j G End(6>5CA) ®M[^/], and also the corresponding 
monodromy matrix 



T A ^a)X M = L Aj *(C/k) . . .L A;/ (C/^)? 2aDA %,/] • 



(21) 



The monodromy matrices (15) and (21) are the basic ingredients in the definition 
of the fermionic operators entering the exponential form of the density matrix. In 
addition to the monodromy matrices we shall need certain operators acting on the 
spin. We introduce the adjoint action of the z-component of the total spin on Af^/j, 



'[*,*] 

and the spin reversal operator defined by 

$ x [k,l] 



(22) 



(23) 



Somewhat loosely we say that an operator X^j has spin s if 



k,l] 



sX\ 



k,l] 



(24) 



We call X^/j spin-reversal even if JX^ = X^/j and spin-reversal odd if JX^/j 



-%/]• 

The annihilation operators b(£, a) and c(£, a) as well as a couple of other interest- 
ing fermionic operators (see [6]) can be derived from a master object k(£, a) defined 
as 

k(£,a)X [M] =tr A;a {o+T a (£^^ (25) 

Here the bosonic trace is understood in such a way that it agrees with the trace in 
the bosonic representation W + = ©£>oC|0) for \q\ < 1 but not a root of unity (with 
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q D \k) = q k \k),a\k) = (1 — q 2k )\k— l),a*\k) = \k+ 1)). For more details see appendix 
A of [6]. 

We shall give a description of the operators b(^, a) and c(^, a) in terms of certain 
residua of k(^, a) which is appropriate for implementing them on a computer. Let us 
fix an operator X[ kJ \ of spin s. Note that \s\ <l — k+\ = n. The following statements 
about the analytic structure of k(^, a) were worked out in [6]. 

(i) t;- a+s - l k(t;,a)X M is rational in t, 2 , 

(ii) Its only singularities in the finite complex plane are poles at C, 2 = ^ 2 ,g ±2 ^ 2 , 

(iii) £- a -* +1 k(C,a)X [M] is regular at t, 2 = oo. 

It follows that k ska i(£,cc)X[ fe;/ ] = ^ a_s ' _1 k(^,a)X[ fc /] is regular at C, 2 = if s < and 
has at most a pole of order s in ^ 2 at C, 2 = if s > 0. Furthermore k s k a i Oi)X^ ~ C, 
for C, — > oo. Hence, 

k s kal(^,CX)X[ fc;/ ] = 



-2 



' « o (e 7a) * 

.7=le=0,± , 9 « Sy y=l 



X 



[k,l] 



(26) 



(where the second sum over j is zero by definition for s < 1). The residua p \ '(a) 
and the Laurent coefficients K,(a) determine all those operators that are needed for 
vanishing magnetic field. They are defined by equation (26). 

First of all there are the operators c(£, a), c(£, a) and f(£, a) defined in section 2.7 
of [6] . In terms of the residua and Laurent coefficients they read 



c(C,a)X [M = ^ +1 £ 



C 2 Ht P 



) Pf («)%,/] 



c(C,a)X [M = ^ +1 H 



^ 2 + y ^(«+-!)pf 



4§ 



f(C,a)Z [M = C 



a+j+l 



■' C~ 2j K i (a)X [t .,i 

„a+i+l-2/' ^-a-5-1+2; 

7=0 *? 

;=le=± £ = 



(27a) 



(27b) 



4* 



where by definition 



pf(«) 



(27c) 



(28) 



j=le=ty± 2^ 2e ^ 2 
In terms of these operators k(£,oc) is decomposed as 

k(£,a)X M = (c(C,a) +c(^,a) +c(^^,a) + f(^,a) -f(rt a))X [M . (29) 
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The operator b(£,cc) that is needed in the construction of the density matrix is 
defined as 

b(C, a) = q- 1 (q^- 1 - q- a+s+l )U(^ -a) J , (30) 

where J is the spin-flip in the adjoint representation. Taking into account that JSJ = 
— 8 we find 

where 

pf («)%/] = 9" 1 (9°-* - <T ") J pf (-a) JX [M , (32) 

for £ = ±. 

For later convenience we introduce the operators 

? -e(a + S + 2)p(e) (a) ^(«+S-2)p(«> (a) 

b »=E > c,-(a)=£ , (33) 

e=± e=± z Sj 

for 7 = 1 , . . . , n. In terms of these operators we have 

b& a)X [M = r a - s+1 £ • b/a)X [M , (34a) 

c& a)X [M = C a+i+1 £ • c,(a)X [M . (34b) 

Let to(/ii,/i2;oc) denote the function that determines the physical part of the prob- 
lem. 0)(jui,jU2;cx) is a known function of the temperature and the magnetic field (see 
section 3 of [2]). Here we only need its property Cfl(//2>j"i;oc) = (Q(jii,fi2', ~ °0- We 
also set o5(jUi,jU2; ct) = co(//i,//2;oc)(£2/£i) a (recall that = e^) which has the same 
symmetry. Then, by definition, 

r dC 2 r dt 2 

^i («)%/] = - J T ^2 j T ^| "O"! ' ^ a)b(Ci , a)c(^ 2 , a - 

= - £ [©(^XyjaJb^aJcyCa-lJ + S^^sa^CaJc/Ca-l)]^ (35) 

l<i<j<n 

(where Xj = ln^y) for all X^/j with S(X[ fc ;]) = 0. The closed contour T in (35) encircles 
all the t, 2 (and excludes the singularities of (0). In the second equation we used that 

b f (a)c/(a-l) = 0, for i = l,...,n. (36) 
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For the physical correlation functions the limit a — > remains to be done. We shall 
write = Qi (0). We claim that 

Q i*M = - E Hh^tj + a '(^j) a 7j] x [k,i], (3V) 

l<i<j<n 



where 



a+ = lim [b,-(a)c 7 -(a- l) + by(a)c,-(a- 1)] , (38a) 



oc^O 



ar. = lim a[b/(a)c y (a- 1) -b 7 -(a)c,-(a- 1)] , (38b) 



and 



<o(Xi,\j) = c«)(^-,\,;0) , co'(^,X 7 -) = a a co(^,X 7 ;a)| a=0 . (39) 

The existence of the limits in (38) was proved in [2]. 

As was already mentioned in the introduction to this section, the operator h(£, a) 
introduced in [2] is a less well understood than b(£, a) and c(£, a). Here we adapt the 
definition of [2] to the novel conventions of [6] and set 

h^a)X M = (l-q a )t^ a {a+T a ^a)T A (^ (40) 

This operator has some similarity with the operator k(£, a) introduced above. The only 
differences are the prefactor ( 1 — q a ) and the insertion of an extra bosonic annihilation 
operator on the right hand side. For this reason its analytic structure as a function of 
£ is rather similar to that of k(£,cc). Following the same reasoning as in section 2.5 
of [6] and using that [§,h(£,cc)] = one sees that for an operator X^ of spin s 

(i) £- a +*h(C,a)X [M is rational in t, 2 , 

(ii) its only singularities in the finite complex plane are poles at C, 2 = l 2 , q ±2 £, 2 j, 

(iii) £- a -*+ 2 h(£, a)X M is regular at t, 2 = <*>. 
It follows that h(£, a)X^ is of the form 



%,/]• (41) 



We denote the function that determines the physical part of ^2 by <PG"; oc) (see [2]) 
and set cp(/j; a) = cp(ju; a)^ a . Then, for every X^ /j of spin 0, 



a 2 X M = - lim / r ^^i;a)h(Ci,a)X [M 



lim £ <p(X y -; a)h ; (a)X M = - £ y(Xj)hjX M , (42) 



10 



where 

h;=limh 7 -(a), q>(Xy)=9(Xy;0). (43) 

Since only the residua at the simple poles at ^ enter the formula for the correlation 
functions, there is a certain degree of arbitrariness in the definition of the operator 
h(C,a). 

Equations (37) and (42) can be used in order to calculate short-distance correlators 
explicitly on a computer. Examples will be shown in section 5 below. The problem 
clearly divides into two separate tasks. The first task is the calculation of Clfj and 
hj. For this purpose one has to generate the operators k(£,cc) and h(£,oc) for a finite 
number of lattice sites and then calculate their residua. The second task is the actual 
calculation of the functions Cfl(//i,//2)> ®'(MiiM2) and <p(fi). They will be described in 
the following section. 

Still, the knowledge of the operators Cli and Q.2 in the form (37) and (42) would not 
be very useful, if these operators would not be nilpotent, which causes the exponential 
series in (9) to terminate after finitely many terms. As was first shown in [3] we have 

V-lP-u = ^tPu = V-i,P-u = , if {/, j} n {*, /} + © , (44a) 

Hp nj/] = Hp a u] = Hp = ° , if {/, j} n {*, / } = (44b) 

which is a consequence of the fermionic nature of the operators b(£, a) and c(£, a). It 
follows for chain segments of length n that 

a^ +1 =0, (45) 

where [x] stands for the integer part of x. Again for Q.2 our knowledge is so far re- 
stricted to what be learned from explicit examples. We calculated Clf-, Q.^- and h, 

l iJ l -J J 

explicitly in the inhomogeneous case for n = 1, 2, 3, 4. The reader might wonder why 
we do not present the matrices here. The reason is that they are very big. Already 
for n = 4 they are so big that they can not be handled any longer by an all-purpose 
computer algebra program like Mathematica. We used the program FORM [19] in- 
stead which turned out to be efficient. We found that for n up to 4 the hj mutually 
anticommute 

h^ + h^h y = (46) 
and that they commute with the Q*-, 

Hp h k]=°- ( 47 ) 

We also have the stronger property 

a^h k = h k ^j = if ke{ij}. (48) 

We conjecture that (46)-(48) are true in general. This would mean, in particular, 
that 

[a h a 2 ] = o (49) 
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and that Q.2 is nilpotent, 

ftf = 0. (50) 
As a further consequence the exponential form (9) would factorize, 

e n 1+ £2 2 = e £2 le n 2 = e "i( 1 + ft 2 ) = (i + Q 2 ) e n i. (51) 

We verified this explicitly for n = 1,2,3,4. 

4 The physical part of the construction 

As we have seen above we need the functions (0(^1, ,1/2), Cfl'Gui,/^) and (p(/j) for the 
description of the physical correlation functions in the limit a — > 0. A relatively sim- 
ple description of these functions in terms of an auxiliary function a and a general- 
ized magnetization density G was given in [2]. We call it the a- formulation. The 
a-formulation is useful for deriving multiple integral formulae [9] and for studying the 
high-temperature expansion of the free energy and the correlation functions [18]. For 
the purpose of the actual numerical calculation of thermodynamic properties [15, 16] 
or short-distance correlators from the multiple integral [7] one has to resort to a differ- 
ent formulation we refer to as the bb-formulation. It needs pairs of functions, but the 
defining integral equations in the massless regime |A| < 1 involve only the real axis as 
integration contour and become trivial in the zero temperature limit. It turns out that 
the bb-formulation is numerically extremely stable (see below). 

Switching from the a- to the bb-formulation is a rather standard procedure which 
basically requires the application of the Fourier transformation (see e.g. [7]). For this 
reason we present only the result. We restrict ourselves to the massless regime |A| < 1 
and set 7 = — ir| e K. Let us define a pair of auxiliary functions as the solution of the 
non-linear integral equations 

, r/ n %h 27L/sin(y) f°° dy w f1 fl 

lnb(x) = -— — - - — , , + / -^-F(x-y)ln(l + b(y)) 

w 2(31-7)7' r7ch(7u/7) }-<*, 2% v y ' v yyn 

^F(x-y + TT)m(l+%)), (52a) 
- lit 



f 



, 7-/ n ™ h 27i7sin(7) r 00 dy^. _ 



271 

with kernel 



^F(x-y-Tl-)ln(l + b(y)) (52b) 



r~ sh (a-y)k)e lkx 

F (x)= M 7^ — / -r- . (53) 

J— 2sh((7i-7)|)ch(|) 

Note that the physical parameters temperature T, magnetic field h, and coupling J enter 
only through the driving terms of equations (52) into our formulae. In this formulation 
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(52) is valid for < y < 7t/2 which means < A < 1. In order to access negative A 
we have to reverse the sign of J and conjugate the local operators whose correlation 
functions we are interested in by o z on every second site of the spin chain (see [21]). 

Except for the auxiliary functions b and b we need two more pairs of functions 
gff^ and g'^ in order to define (0 and go'. Both pairs satisfy linear integral equations 
involving b and b, 

gi +) (*) = fsech(^) 

[°° dy F(x-y) (+), , f°° dy F(x- y + r\~) r_), N /e . . 

♦/..HITF^ M ~Lrn 1+r . M 'W. <*■> 

g<-»W = f S ech(lfcii) 

p dy F(x-y) H f- dy Fjj-y-TT) (+ , 

+ /..srrr^«' w-Ls-i+Fi-* (>,) ' ' 



and 



g'i +) (x)=(f(x- / i)-f)sech(fe)) 



gl } (*)= (f (*-,/) + §)sech(^ 



sin(fcc)sh(£ )ch((5-y)Jfc) 
D(x) = / dfc — — K2 > - KK \ [' ' . (56) 
1 ' J— 4sh 2 ((7i-y)f)ch 2 (f) 

The functions Cfl(//i,//2), Co'Gui,^) and (p(//) that determine the explicit form of 
the inhomogeneous correlation functions of the XXZ chain can be written as integrals 

involving b, b, gj^ and g'^. The function 



where 



<pG«) = /_ 



cbc 



-co 2(7i — y)i 



1 + b !(*) 1 + b 



(57) 
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also determines the magnetization 

m(r,fc) = -i(p(0) (58) 
which is the only independent one-point function of the XXZ chain. The function 

sh((ji-y)|) cos(&(//i-// 2 )) 



-« 

j — < 



dk- 



^sech ^ 



ish(f)ch(¥) 

„(+) 



iS'/'i ( x ) 



■ + 



i+b-^x) i+r\x) 



with 



K(/u)=cth(/j-r\) -cth(// + Ti) 
also determines the internal energy 

e (r,/i)=7(^_ 1 a) + a^ 1 ^ + A(o^ 1 o}-l)) r ,, = 7(sh(ri)a)(0,0)-A) 
of the XXZ chain. The function 0/(^1,^) is defined as 

Ysin(A;(/ii -^2)) 



(59) 

(60) 
(61) 



(o / ( w , A / 2 ) = ^(+)( Ail - jU2 ) + 



dfc 



2ith(f)ch^) 



+ £^-, 2 )sech(^) 



l + b l + b 



(62) 



where 



^(+)(^) = cthOi-ti) + cth(^ + Ti) , /^(x) =4 ±) W=F^i ±) W. (63) 

For the function a/(0, 0) itself we did not find a simple interpretation in terms of ther- 
modynamic quantities. Its derivative Qd f x — 3 jU (0 / (^, 0)1^=0, however, appears in the ho- 
mogeneous density matrix Z) 2 (r,/z) and, hence, is related to the neighbour correlators 
(see [2]). 

From the symmetry Cfl(//i,// 2 ,a) = Cfl(// 2 ,//i, —a) (see [2]) we conclude that 



0>G"l,/*2) = ®(M2,Pl) , (if (MUM) = -®'(M2,Pl) 



(64) 



This property was crucial for the derivation of (37). 

For the calculation of the physical correlation functions one has to perform the ho- 
mogeneous limit. This limit is rather singular, because the coefficients coming from the 
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algebraic part in general have poles when two inhomogeneities coincide. These poles 
are canceled by zeros stemming from certain symmetric combinations of the functions 
(o(A,i,A,2), Cfl'(A,i, A2) and (p(A) in the numerators. In order to perform the limit one 
has to apply 1' Hospital's rule, finally leading to polynomials in the functions (o(0,0), 
g/(0,0), (p(0) and derivatives of these functions with respect to the inhomogeneity 
parameters evaluated at zero. We shall denote derivatives with respect to the first argu- 
ment by subscripts x and derivatives with respect to the second argument by subscripts 
y and leave out zero arguments for simplicity. Then e.g. G>' xyy = d x dy(o'(x,y) \ x ,y=o etc. 

For the examples in the next section the non-linear integral equations for b and b 
as well as their linear counterparts for and g'^ were solved iteratively in Fourier 
space, utilizing frequently the fast Fourier transformation algorithm. The derivatives 
of g^ and g'^ with respect to /j, needed in the computation of the respective deriva- 
tives of (0, G)' and (p, satisfy linear integral equations as well, which were obtained 

as derivatives of the equations for g^ and g'^\ Taking into account derivatives is 
particularly simple in Fourier space. 



5 Examples of short-distance correlators 

Here we mostly concentrate on the longitudinal and transversal two-point functions 
{<5\<5 z n )T,h, (of G^Tji for n = 2,3,4. The algebraic partforn = 2,3 was already calcu- 
lated in the more general inhomogeneous case in our previous paper [2]. 
In the inhomogeneous case for n = 2 we obtained 

= - ^ tri2 ( [co( Ai , A 2 ) ^ 2 + ( Ai , A 2 ) a^) 

= cth(r [ )(D(XiA2) + Cth( ^~ X2) a) / (X 1 A2) (65) 



and in the homogeneous limit Xj — > 0, 



{c\o z 2 ) T ,h = cth(Ti)co+^. (66) 

Note that the limit exists, because Cfl'(Ai,A2) is antisymmetric. Similarly for the 
transversal correlation functions, 

itrizfe" 1 (1+^2)0102) = ^tri 2 (Giofo9 

= - i tri 2 ( M) &t,2 + &/ M 

ch(Ai-A 2 ) » v ch(ii) , 

-C0(A,i,A,2)- — r— (0(Ai,A 2 ) (67) 



2sh(ri) v 2rish(Ai-A 2 ) 
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which in the homogenous limit becomes 

( ^^-"2lhM — 2iT- m 

The two-point functions (o\o z n )T,h> {<s\<5 x n )T,h are even under spin reversal. There- 
fore we do not expect any non-vanishing contribution from ^2, which would come 
with functions 9, i.e. odd functions of the magnetic field h. This is indeed what we ob- 
serve from our calculation for n = 1, 2, 3,4. The above two-point correlation functions 
depend only on (0, a/ and their derivatives. As an example of a correlation function 
which is of no definite symmetry under spin reversal and which does depend on cp and 
its derivatives we shall discuss the emptiness formation probability below. 

We obtained the formula for n = 3,4 in the same way as for n = 2 by first calculat- 
ing the inhomogeneous generalization of the correlation function and performing the 
homogeneous limit at the end. The formula for in the inhomogeneous case 

was presented in [2] . For space limitations we do not repeat it here and also do not 
show the inhomogeneous formulae in the remaining cases. Instead we merely list the 
results for the physical correlation functions. For n = 3we already obtained in [2] 

/ z z\ o .urn ^ , , th(ri) ((0^-2(0^,) sh 2 (r|)a4 

(°i03>r,A = 2cth(2r|)a)+ + 4 , (69a) 



1 m ch ( 2r l) f ,y 



x 



ch(2ri)th(ri)((o xx -2(o xy ) sh 2 (ri)a)^ 

8 + ~8^ • (69b) 

The length of the formulae grows rapidly with the number of lattice sites. For n = 4 
we found 

(°l°4>r,A 1 



768<7 4 (-l+<7 6 )(l+<7 2 )r| 2 

1 3 8 V ( l + ^ 2 ) 2 (5 - 4^ 2 + 5^ 4 ) ri 2 co 

- 8 (l + q 4 (52 + 64q Z - 234q 4 + 64q 6 + 52q S + q l2 ^j ) T| 2 (0, 

+ 192<7 4 (-1 + q 2 ) 2 (1 + 4q 2 + q 4 ) V[ 2 (O yy 
+ (- 1 + q 2 ) 4 ( 1 + q 4 ) ( 1 + 4^ 2 + q 4 ) T| 2 [-4(0 W + 6(0 xxyy j 

- 16Sq 4 (-l-q 2 + q 6 + q S ^j r\(ii' y 
+ 16 (- 1 + q 2 ) 3 (l + 6q 2 + 1 \q 4 + 1 \q 6 + 6q 8 + <? 10 ) 1^ 
- 2 (- 1 + q 2 ) 5 (l + 2q 2 + 2q 4 + q^ T|Co'. 



' xxyyy 
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+ 8 (- 1 + q 2 ) 3 ( 1 + q 2 ) ( 1 + 6q 2 + 34q 4 + 6q 6 + q^r[ 
+ (-1 - Aq 2 - 22q 4 -12q 6 + 12q 10 + 22q 12 + Aq u + q 16 ^j T| 



-6(0 



yy 



+ 12(%C0x>, + 4co v (0 



•yvjyyy 



12(0y(0 



'xyy ' 



AdXHxyyy + 6(0(0. 



'xxyy 



+ \6(-\+q 2 ) 4 (\+q 2 ) 2 (l+q 2 + q 4 )T} 



QOyyQO'y ~ QOyQO' yy + GKti' X yy 



+ ( - 1 + q 4 ) 2 ( 1 + 5q 2 + 6q 4 + 5q 6 + q 8 ) T| 4(£> xyyy G)'y - 6(£> xxyy &'y - 2co yxy Cfl / 

+ 6Cdxyy(i)' yy + 2(Hyy(Sl yyy - 4(O xy Q0' ' yyy - 6(DyyQd' x yy + 4 QOy QO' x yyy ~ 2(0^ ' XX yyy 



+ 3(-l+q 4 ) 3 (l+q 2 + q 4 ) 



(0 yyy (H x yy (0 yy& xyyy (0 y(0 xxyyy 



and in the transversal case, 



(70) 



T,h = 



3072q 5 (-l+q 6 )r\ 2 
|-768^ 6 (l + 10^ 2 + ^ 4 )ri 2 (o 

+ I64 2 (- 1 + q 2 ) 2 (3 1 + 56^ 2 - 30q 4 + 56q 6 + 3 lq s ^j T| 2 (0 
- 96q 2 (-l+q 2 ) 2 (3 + 5q 2 - Aq 4 + 5q 6 + 3<? 8 ) y\ 2 ® yy 



+ q 2 (-l+q 2 ) 4 (l+4q 2 + q 4 )TY 

+ \92q 2 (-3-q 2 -q 4 + q*+ q 10 + 3q 12 ) T\0)' y 
+ 8 (- 1 + q 2 ) 3 (l - 12^ 2 - 25q 4 - 25q 6 -12q s + q 10 ^j Tf]0i' 

+ 2 (- 1 + q 2 ) 5 (l + 2<? 2 + 2q 4 + /) T\d xxyyy 

+ \6q 2 {-l+q 2 ) 3 (\l + lq 2 + lq 4 + 17<? 6 ) T] 2 

+ q 2 (-5 - 4<7 2 -I3q 4 + \3q 8 + 4q w + 5q 12 ) T] 2 



8 ($xyyy 1 2 (Oxxyy 



(0(0 xy - QOy 



120)^ y -24(Oyy(0 X y 



S&y&yyy + 24(tiyG) x yy + SCOCO^yy — 1 2&& XX yy 



+ 8 (- 1 + q 2 ) 4 ( 1 - 9q 2 - %q 4 - 9q 6 + <? 8 ) T| CQyyCfl'j, - CQyCfl'^ + (D(D f xyy 
+ (- 1 + q 4 ) 2 ( 1 + 5q 2 + 6q 4 + 5q 6 + <? 8 ) T| -4(0^(0^ + 6(0^(0^ 
+ 2(0^(0'^ — 6cQ; W CQ / xy — 2(o xy (0 / ;TO , +4cQxyCQ / xxy + 6co xy co'x Xy — 4co y Cfl'x Xxy + 2&d xxyy y 



3 

+ 3 (- 1 + q 4 ) ( 1 + q 2 + q 4 ) -Cfl'^Cfl'^ + Co'^G)'^ - (0^(0^^ 



• (71) 




Figure 1 : Two-point correlators for n = 2, 3, 4 as a function of temperature for various 
values of A in the massless regime at zero magnetic field. 
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.464) 
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0.588818 
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(0, 


.355) 


-0.7 


0.412795 


(0.413) 


0.316321 


(0, 


,318) 


0.262606 


(0, 


.264) 


-0.8 


0.262355 


(0.265) 


0.211402 


(0, 


,215) 


0.179803 


(0. 


.184) 


-0.9 


0.129195 


(0.137) 


0.111127 


(0. 


.118) 


0.098055 


(0. 


.104) 



Table 1: Numerical values of the cross-over temperature 7o(n;A) on the infinite chain 
and for a chain of 18 sites in units of J/2 (values in brackets, taken from [8]). 

"The prefactor on the rhs of equation (70) has a pole as a function of q when q is a third root of unity 
corresponding to A = —0.5. We estimated the value of 7b(n; A) for A = —0.5 by interpolating between 
two values closely above and below —0.5. 

Using the representations of the previous section for the functions (0 and a/ we can 
determine high-precision numerical values for the various two-point correlators. Fig- 
ures 1-3 show selected examples. For n = 2, 3, 4 the longitudinal correlation functions 
(a\o z n )T,h are exhibited in the left panels, while the right panels show the transversal 
correlation functions (o^o^tm- We observe a rich scenario for their temperature and 
field dependence, in particular non-monotonous behaviour at intermediate tempera- 
tures. The numerical precision in all figures is between 3 and 4 significant digits. 

In figure 1 we display our data for vanishing magnetic field and various values of 
A G ( — 1,1). The transversal correlations are alternating in sign with n for all values 
of A. But for the longitudinal correlations their qualitative behaviour depends on the 
sign of A. They alternate for positive A. For negative A they are always negative at 
sufficiently low temperatures, change their sign at a 'crossover temperature' 7o(n;A), 
reach a positive maximum and tend to zero from above in the high temperature limit. 
This phenomenon was studied numerically on finite chains of 16 and 18 sites in [8] and 
was interpreted there as a 'quantum-to-classical crossover'. In [8] it was estimated that 
the crossover temperature Tq{h; A) remains finite in the thermodynamic limit. From our 
formulae we can calculate it with high precision (see table 1). 

Note the peculiar low temperature behaviour close to A = — 1, where one would 
expect simultaneously isotropy ((g^o^) = (— 1)" +1 (c\c*)) as well as full polarization 
((GjO^) +2(— l) n+l (c\o*) = 1). However, this is not seen, not even for A as close to 
the isotropic point as —0.995, since the full polarization is still forced into the xy-plane 
by the slightly dominant transversal exchange. At elevated but not too high tempera- 
ture the system is still close to full polarization, but the small anisotropy is irrelevant 
and longitudinal as well as transversal correlation attain values close to ±1/3. At much 
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Figure 2: Two-point correlators for n = 2, 3, 4 and A = 1 / v2 as a function of temper- 
ature for various magnetic fields. 
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higher temperatures all correlations tend to zero. As was explained in [8] the negative 
sign of the longitudinal correlations at low temperatures is a pure quantum effect. 

In figure 2 we plotted the correlators for a fixed value of A = 1 / \fl for various mag- 
netic fields. The curves show non-monotonous behaviour for intermediate fields. In 
particular, a sign change is possible for positive A if the magnetic field is switched on, 
as can be seen in the upper and lower left panels. In general, we observe a competition 
between the effect of the magnetic field, which tends to align the spins in z-direction, 
the longitudinal exchange interaction in the Hamiltonian, which enforces antiparallel 
alignment of neighbouring spins, and the transversal exchange interaction, which can 
be interpreted in terms of quantum fluctuations. The temperature sets the scale for 
the relative and absolute relevance of these terms. At low temperatures the quantum 
fluctuations dominate, decreasing e.g. the value of (<3\o z 3 ) from 1 to a smaller positive 
value. When increasing the temperature, at intermediate field the quantum fluctuations 
are suppressed and the field is less hindered to align the spins, before at very high 
temperature the thermal fluctuations destroy all order. 

For large enough fields (light blue curves) the correlation functions reach satura- 
tion. This can be also observed in figure 3, where the correlators are shown as functions 
of the magnetic field for relatively low temperature, T/J = 0.1, and where the satura- 
tion roughly sets in at its zero temperature value h = (1 + A)47. An interesting feature 
is the decrease of (crjOg) for small fields which we interpret as due to a weakening of 
the antiferromagnetic neighbour correlations. 

As an independent test of our results, we numerically calculated the correlation 
functions for finite systems. Since the Hilbert space dimension grows exponentially 
with the chain length, this required us to characterise the spectrum and eigenvectors 
of high-dimensional sparse matrices. The complete diagonalisation of such matrices 
is too time consuming or simply impossible due to memory limitations. We therefore 
used Krylov space methods, in particular, Chebyshev expansion [20]. 

Given the D eigenstates \k) of the system, we split the expectation value of an 
operator A into a low-energy and a high-energy part, 



Z = Tr(e-"/ r ) = £ e- £ ^ r + £ e^/ r , (72) 



<A> = i T r(Ae-"/ r ) = I 



C-l D-l 
£=0 k=C 

D-l 

" ' ' (73) 



C-l D-l 

£ (k\A\k) e^/ r + £ (k\A\k) e-^/ r 

k=0 k=C 



where H is the Hamiltonian of the system, Z the partition function, and the energy 
of the eigenstate \k). The C low-energy states and the corresponding matrix elements 
can be calculated exactly with the Lanczos algorithm [17]. The contribution of the 
remaining states is estimated via Chebyshev expansions of the density of states, 



£ 1 e -V^ = f p( E )c- E / T dE, (74) 
k=c J 




Figure 3: Two-point correlators for n = 2, 3, 4 and T/J = 0. 1 as a function of magnetic 
field for various values of A. 
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Figure 4: Two-point correlators for n = 4 and A = 1 / \/2 as a function of temperature, 
comparison with data from direct diagonalization of finite chains. 



D-l 



1 



p(E)=^5(E-E k ) 

fe=c 7iv* z -(^-^r 

/ig=g?Tr(Pr B [(H-£)A]) with P = 
and of the diagonal matrix elements of A, 

D-l 

E 

k=C 



: i-Ei*x*i, 

k=0 



E-E - 
s ' 



(75) 
(76) 



Y,{k\A\k)z- Ek ' T = [ a(E) eT E l T dE 

D-l 

a(E)= £{k\A\k)8(E-E k ) 



1 



k=c ny/£-(E-E)* 
V a n =g%Tr(PAT n [(H-E)/s]). 



n=l 



E-E-- 



(77) 

(78) 
(79) 



Here T n are the Chebyshev polynomials of the first kind and E, s denote scaling factors 
that shift the spectrum of H into the domain of the polynomials [—1,1]. The coef- 
ficients account for a special kernel [12,20] which improves the convergence of 
the truncated Chebyshev series. The traces defining the moments /uf, and /u" can be 
estimated stochastically, 

R-l 

Tr(X) « 2jf ( r \ x \ r ) with R < D random states |r) . (80) 

r=0 

Thus the most time consuming part of the moment evaluation is the recursive calcula- 
tion of T n [(H — E)/s] | r) for a few dozen states | r) . The whole procedure is feasible for 
Hilbert space dimensions of 10 7 to 10 8 . Employing S z and momentum conservation, 




Figure 5: Emptiness formation probability for n = 2, 3 and A = 1/ y/2 as a function of 
temperature for various magnetic fields. 



with moderate effort we can reach a chain length of L = 32 (about 430 CPU hours on 
a 2.66 GHz Xeon Cluster). 

We also performed tests against other exact results. We compared the high temper- 
ature expansions of equations (70), (71) with the high temperature expansions for the 
multiple integrals from [10] up to third order in l/T. We compared with the known 
ground state results of [13, 14]. We compared the zz-correlation functions in the XX 
limit with the explicit result (see e.g. [1 1]), and we computed the correlation functions 
for A very close to 1 and compared with the results obtained for the isotropic limit and 
w = 2,3in[l]. 

As an example of a correlation function which has no definite symmetry with re- 
spect to the reversal of all spins and hence must depend on cp and its derivatives we 
consider the emptiness formation probability P(n) = 2~"(n" = i(l +<^))r,/i- We know 
from [2] that 

1 <P cth(ri)(o co' 
P(2) = --| + ^^ + ^ (81) 

and 

P n) = I _ + ( 1 +^ 2 + ^ 4 ) (0 , (-l+<? 2 ) jj%-2tQgj _ 3ro; 
U 8 8 2(-l+4 4 ) 32(1+^ 2 ) 8r| 

(-1+4 2 ) <yy 3 (1 + (-WCpXt + 2(Oy(p X + (Oyy^ ~ 2(0^) 



128^ 32(-l+q 



2\ 



(l + l0q Z + q 4 ) ((0^9-0)^ + 0)^) 
+ 128^ ■ (82) 



P{2) and P(3) are displayed in figure 5. They show the characteristic high temperature 
behaviour P(n) ~ 2~" and saturate for small temperature and h > (1 + A) 47. Like the 
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two-point functions they exhibit non-monotonous behaviour at intermediate magnetic 
fields. 

6 Conclusion 

The aim of this work was to demonstrate that short-distance correlation functions of 
the XXZ chain can be computed efficiently from the exponential form of the density 
matrix proposed in [2]. We found that an accurate computation is possible for arbitrary 
temperatures and magnetic fields in the massless regime —1 < A < 1, even beyond 
the phase boundary to the fully polarized state. This has to be contrasted with our 
experience [7] in computing the same correlation functions from the multiple integral 
representations obtained in [9, 10], which appeared to be much harder. In fact, in [7] we 
were unable to proceed beyond n = 3. In the present approach the difficulty for larger n 
comes mainly from the algebraic part of the calculation. We believe that ranges of the 
order of n = 10 may be reached if one works directly with the homogeneous transfer 
matrices (which would mean to calculate higher order residua in (9)). 
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